#include "Polynomicon.h"
#include "FourierTron.h"
#include <iostream>
#include <cmath>


using namespace std;
using namespace SeamStress;


void Polynomicon::modPoly(const refloat_vec& input, const refloat_vec& divisor, refloat_vec& remainder)
{
  if(input.size() < divisor.size())
  {
    remainder = input;
    return;
  }
  
  refloat_vec rev_b;
  for(long int i=(divisor.size()-1);i>=0;i-=1)
  {
    rev_b.push_back(divisor[i]);
  }
  refloat_vec rev_b_inv;
  this->reciprocalPowerSeries(1+input.size()-divisor.size(), rev_b, rev_b_inv);
  modPoly(input, divisor, rev_b_inv, remainder);
}


void Polynomicon::modPoly(const refloat_vec& input, const refloat_vec& divisor, const refloat_vec& rev_b_inv, refloat_vec& remainder)
{
  if(input.size() < divisor.size())
  {
    remainder = input;
    return;
  }
  
  refloat_vec rev_a;
  for(long int i=(input.size()-1);i>=0;i-=1)
  {
    rev_a.push_back(input[i]);
  }
  refloat_vec rev_q;
  this->multiplyPoly(rev_a, rev_b_inv, rev_q);
  rev_q.resize(1+input.size()-divisor.size(), retempval1);
  refloat_vec q;
  for(long int i=(rev_q.size()-1);i>=0;i-=1)
  {
    q.push_back(rev_q[i]);
  }
  this->multiplyPoly(q, divisor, remainder);
  retempval1.set(0);
  remainder.resize(divisor.size()-1, retempval1);
  retempvecpointer1 = &remainder;
  retempvecpointer2 = &input;
  pins.sewStraight(&Polynomicon::modPolyThread1, nthreads);
}


void Polynomicon::modPolyThread1(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  for(unsigned long int i=w, loopsizei=retempvecpointer1->size();i<loopsizei;i+=nthreads)
  {
    (*retempvecpointer1)[i] = (*retempvecpointer2)[i] - (*retempvecpointer1)[i];
  }
}


// R = Q(a*x + b*x^2)
void Polynomicon::composeQuadraticPoly(const refloat_vec& Q, const refloat& a, const refloat& b, refloat_vec& R)
{
  // z == ax + bx^2 = b*( (a/b)x + x^2 ) == b*( c*x + x^2 )
  // k == z/b .  z = k*b
  
  refloat_vec cvec1 = Q;
  this->scalePoly(cvec1, b);
  // cvec1 is a polynomial in k
  // k = ( c*x + x^2 ) = (x + c/2)^2 - c^2/4
  // y == ( x + c/2 )^2
  // k = y - c^2/4 == y - y0
  
  refloat c(0, prec);c = a/b;
  refloat y0(0,prec);y0 = -c*c;y0.mul_power_2(-2);
  this->shiftPoly(cvec1, y0);
  
  // cvec1 is a polynomial in y
  // w == ( x + c/2 )
  // y = w^2
  
  refloat temp1(0,prec);
  refloat_vec cvec2;cvec2.assign(2 * (cvec1.size()), temp1 );
  for(unsigned long int i=0, loopsizei=cvec1.size();i<loopsizei;i+=1)
  {
    cvec2[2*i] = cvec1[i];
  }
  
  // now cvec2 is a polynomial in w
  
  // x0 == -c/2
  // w = x - x0
  
  refloat x0(0,prec);x0 = c;x0.mul_power_2(-1);
  this->shiftPoly(cvec2, x0);
  
  // now cvec2 is a polynomial in x
  
  R = cvec2;
  R.resize(Q.size(), temp1);
}


// R = Q(P)
void Polynomicon::composePolyRecursive(unsigned long int n, const refloat_vec& Q, const refloat_vec& P, refloat_vec& R)
{
  R.clear();
  if( Q.size() == 2)
  {
    R.resize(P.size(), P[0]);
    for(unsigned long int i=0;i<P.size();++i)
    {
      R[i].set(P[i]);
    }
    vectorMultiply(R, Q[1]);
    R[0].add(Q[0]);
    retempval1.set(0);
    if(R.size() > (n+1))
    {
      retempval1.set(0);
      R.resize((n+1), retempval1);
    }
    return;
  }
  if(Q.size() == 1)
  {
    R.push_back(Q[0]);
    retempval1.set(0);
    if(R.size() > (n+1))
    {
      retempval1.set(0);
      R.resize((n+1), retempval1);
    }
    return;
  }
  
  unsigned long int qn = Q.size();
  refloat_vec Qtemp;
  Qtemp.resize(qn,Q[0]);
  for(unsigned long int i=0;i<qn;++i){Qtemp[i].set(Q[i]);}
  unsigned long int qn2 = 1;
  while((qn2+1) < qn){qn2<<=1;}
  qn2+=1;
  retempval1.set(0);
  if(qn2 > qn){Qtemp.resize(qn2,retempval1);}
  
  refloat_vec Q1;
  Q1.resize((qn2-1)/2, Q[0]);
  for(unsigned long int i=0;i<Q1.size();++i){Q1[i].set(Qtemp[i]);}
  refloat_vec R1;
  composePolyRecursive(n, Q1, P, R1);
  if(R1.size() > (n+1)){R1.resize(n+1,retempval1);}
  
  refloat_vec Ptemp1, Ptemp2, Ppower;
  Ppower.resize(P.size(), P[0]);
  for(unsigned long int i=0;i<P.size();++i){Ppower[i].set(P[i]);}
  unsigned long int power = 1;
  while( power < ((qn2-1)/2) )
  {
    Ptemp1.resize(Ppower.size(), Ppower[0]);
    for(unsigned long int i=0;i<Ppower.size();++i){Ptemp1[i].set(Ppower[i]);}
    Ptemp2.resize(Ppower.size(), Ppower[0]);
    for(unsigned long int i=0;i<Ppower.size();++i){Ptemp2[i].set(Ppower[i]);}
    multiplyPoly(Ptemp1, Ptemp2, Ppower);
    if(Ppower.size() > (n+1)){Ppower.resize(n+1, retempval1);}
    power*=2;
  }
  
  refloat_vec Q2;
  Q2.resize((qn2-1)/2 + 1, Q[0]);
  for(unsigned long int i=0;i<Q2.size();++i){Q2[i].set(Qtemp[i + (qn2-1)/2]);}
  refloat_vec R2;
  composePolyRecursive(n, Q2, P, R2);
  
  multiplyPoly(R2, Ppower, R);
  if(R.size() > (n+1)){R.resize(n+1,retempval1);}
  addPoly(R1, R);
  if(R.size() < (n+1))
  {
    retempval1.set(0);
    R.resize((n+1), retempval1);
  }
}


void Polynomicon::composePoly(const refloat_vec& Q, const refloat_vec& P, refloat_vec& R)
{
  unsigned long int n = P.size() - 1;
  unsigned long int m = (unsigned long int)(sqrt( ((double)n)/log((double)n)));
  if(m == 0){m=1;}
  unsigned long int l = n/m + 1;
  calculateFactorials(l+1);
  refloat_vec Ph;Ph.assign(m+1, P[0]);
  for(unsigned long int i=0;i<(m+1);++i){Ph[i].set(P[i]);}
  retempval1.set(0);
  refloat_vec Pr;Pr.assign(n+1, retempval1);
  for(unsigned long int i=(m+1);i<(n+1);++i){Pr[i].set(P[i]);}
  
  refloat_vec Pr1, Pr2;
  refloat_vec* Pr_cur = &Pr1;
  refloat_vec* Pr_pre = &Pr2;
  Pr1.assign(n+1, Pr[0]);
  for(unsigned long int i=0;i<(n+1);++i){Pr1[i].set(Pr[i]);}
  
  refloat_vec QPh1, QPh2;
  refloat_vec* QPh_cur = &QPh1;
  refloat_vec* QPh_pre = &QPh2;
  composePolyRecursive(n, Q, Ph, *QPh_cur);
  R.clear();
  retempval1.set(0);
  R.assign(n+1, retempval1);
  addPoly(*QPh_cur, R);
  refloat_vec derivative_series;
  differentiatePoly(Ph, derivative_series);
  retempval1.set(0);derivative_series.resize(n+1, retempval1);
  refloat_vec invPprime;
  reciprocalPowerSeries(n+1, derivative_series, invPprime);
  refloat temp1(0,prec);
  refloat_vec Rtemp;
  unsigned long int j=1;
  while(true)
  {
    swap(QPh_cur, QPh_pre);
    differentiatePoly(*QPh_pre, derivative_series);
    multiplyPoly(derivative_series, invPprime, *QPh_cur);
    if(QPh_cur->size() > (n+1)){QPh_cur->resize(n+1, retempval1);}
    multiplyPoly(*QPh_cur, *Pr_cur, Rtemp);
    if(Rtemp.size() > (n+1)){Rtemp.resize(n+1, retempval1);}
    factorialRecipVal(j, temp1);
    vectorMultiply(Rtemp, temp1);
    addPoly(Rtemp, R);
    if(j==l){break;}
    multiplyPoly(*Pr_cur, Pr, *Pr_pre);
    swap(Pr_cur, Pr_pre);
    if(Pr_cur->size() > (n+1)){Pr_cur->resize(n+1, retempval1);}
    j+=1;
  }
}


// sum_{i=0}^{n} a_i*x^i -> sum_{i=0}^{n-1} a_(i+1) * (i+1) * x^i
void Polynomicon::differentiatePoly(const refloat_vec& input, refloat_vec& output)
{
  retempvecpointer2 = &input;
  output.clear();
  retempval1.set(0);
  output.resize(input.size()-1, retempval1);
  thread_var1 = input.size() - 1;
  thread_tot = nthreads;
  if(thread_tot > thread_var1){thread_tot = thread_var1;}
  retempvecpointer1 = &output;
  pins.sewStraight(&Polynomicon::differentiatePolyThreadReal, thread_tot);
}


void Polynomicon::differentiatePolyThreadReal(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  
  for(unsigned long int i=w;i<thread_var1;i+=thread_tot)
  {
    rethread_vec1[w].set(i+1);
    mul( (*retempvecpointer2)[i+1], rethread_vec1[w], (*retempvecpointer1)[i] );
  }
}


